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Abstract 

Monte Carlo simulations based on a Spin-1 Ising Model for binary alloys have been used 
to investigate the non-equilibrium partition coefficient ( k neq ) as a function of solid-liquid 
interface velocity and orientation. In simulations of Si with a second component k neq 
is greater in the [111] direction than the [100] direction in agreement with experimental 
results reported by Aziz et al. The simulated partition coefficient scales with the square 
of the step velocity divided by the diffusion coefficient of the secondary component in the 
liquid. 

Introduction 

During rapid solidification of a binary alloy, the solid-liquid interface is not at equi- 
librium. The partition (or segregation) coefficient therefore differs from that determined 
from the equilibrium phase diagram. Rapid solidification leads to the development of 
novel microstructures and compositions not accessible by standard solidification tech- 
niques. Rapid solidification may be effected by such techniques as laser melting, splat 
cooling or melt atomization. [1, 2] Both analytical and computational models have been 
developed for k neq . The analytical model most extensively compared with experiment is 
due to Aziz. [3, 4] Jackson has investigated this phenomenon using a computational Spin- 1 
Ising model. [5, 6, 7] In addition, orientational dependence of the partition coefficient of 
Bi in Si during rapid solidification was experimentally investigated by Aziz et al. [8] 


Simulations of Solidification of Doped Silicon 

The probability of an atom going from a liquid to solid is given by: [9] 

P t LS = P 0 {K . + h i ) exp ex P H j ( 1 ) 
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is the bond energy between an atom i and an atom j. The energy of the bond between 
two atoms depends only on the types of atoms involved, and the relevant phase diagram. 
(In this paper the liquid was assumed to be an ideal solution, while the solid was assumed 
to be a regular solution.) The possible types of atoms are solid of species a , solid of 
species b, liquid of species a, and liquid of species b, which are represented by "A ,‘B ,'a , 
and ‘b’ respectively. AS is the entropy of fusion, S t , is one when atom i is of type t, k B is 
Boltzmann’s constant and T is the temperature. P Q is a constant that will be determined 
if all the probabilities are to sum to one. Alternatively, if one considers the probabilities 
to be the likelihood of an event in a particular time interval, P 0 will determine the time 
scale. The probability of an atom going from crystal to liquid is given by: 


Pf L = P 0 {t>Ai + <Wexp 


£ 


<6 

2_o 

k B T 


(2) 


Diffusion is represented by an exchange between a liquid atom of type k b’ with a liquid 
atom of type ‘a’ with the following probability: 


P? = PoTSbi 
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Where V is the diffusive jump frequency, Z is the number of nearest neighbors. This will 
lead to a diffusion coefficient given by [10] 
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where a is the jump distance to the nearest neighbor and d is the dimensionality of the 
lattice. 


Partition Coefficient Determination 


Under non-equilibrium conditions the gradient of concentration in the liquid near the 
interface can be very large, resulting in a concentration in the liquid that varies in atomic 
dimensions. Therefore care must be used in the definition of the segregation coefficient. 
A definition of k neq can be based on the solution of the 1-d moving boundary problem: 
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where the concentration (represented by u) of the solid layer is k neq times the concentration 
of the neighboring liquid layer. The equation is differenced where Ax does not approach 
zero, but is determined by the atomic lattice spacing. For silicon the (111) spacing of 
3.13 A was chosen. The resulting difference equation is: 
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Figure 1: Ising model simulation results for silicon at different temperatures, (a) growth rate for the 
[100] and [111] growth directions, (b) Step density in the [111] and [100] directions. 


The subscripts 1, m , and r refer to left middle and right, i is the node number which is 
is increasing to the right, and n is the time interval [11]. The fluxes can be considered to 
be into and out of individual atomic planes. The coefficients to use are given by Table 1 
where B stands for the boundary layer, L a liquid Layer, and S the Solid Layer. B(L) 
represents the boundary with a liquid layer to its right, (L)L represents a liquid layer with 
a liquid layer to its left. 
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Using this procedure a non-equilibrium partition coefficient can be fit to the Monte Carlo 
simulation results ( using the appropriate diffusion coefficient and growth velocity). Fig- 
ure 2 shows a comparison of finite difference calculations and Monte Carlo simulations. 


Scaling of Monte Carlo Simulation Parameters 

The paramaters in the Spin-1 Ising Model are the nearest neighbor bond energies, the 
latent heat of fusion at the melting temperature, the diffusive jump frequency T, and the 
crystal structure. For the present work a diamond cubic structure with a lattice spacing 
of 5.43 A was used. In order for the Ising model to give the proper surface roughening 
temperatures, an entropy of fusion greater than the experimental value must be used. 
The value for A S/k B used in the simulations was 6. Interpretation of the results depends 


3 






Figure 2: The initial concentration profile was the same for all simulations. Concentration profiles of 
the solid were all convoluted with a Gaussian before comparison (convoluted profiles shown). Both series 
of MC simulations were in the [111] direction with ft = 10 and a temperature of 1450 K. The parameters 
for the finite difference run were D — 2xl0 — a velocity of 14.4 tn/s, k neq value of 0.42, and a grid 
spacing of 3.13 A. In this case only the value of k ntq was used as a free parameter for fitting the finite 
difference run to the Ising model simulation results. 


on P 0 in equations 1,2 and 3. By matching the growth velocity in the -simulations to 
experimental growth velocities it was found the time scale could be changed to seconds 
by setting P Q «s 9xl0 13 . ( A single value for P Q was found by fitting experimental data at 
1550 K. In general P 0 should account for a temperature dependent arrival rate at steps.) 
T is defined to be: 

r = QR a (7) 

where R a is the arrival rate at a step, and ft is a parameter which can be adjusted to 
match diffusion data. For data chosen equation 1 simplifies to: 

R a = P 0 exp(—AS/k B ) (8) 


Using equation 4 for diamond cubic and a value of 10 for ft, which is reasonable on the 
basis of atomic considerations [12], we have: 
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Results and Discussion 

Figure 1 (a) shows simulation results for growth of silicon in the [111] and [100] directions. 
The line for the [1 1 1] direction should come in vertically to the melting point of silicon since 
the surface is smooth. Figure 1 (b) shows the step density as a function of undercooling 
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Figure 3; Ising model simulations for the diamond cubic structure in the [111] and [100] directions with 
ft = 10. A horizontal line is drawn at k eq = .001. (a) k neq as a function of normalized velocity, (b) k n , q 
as a function of normalized step velocity. 


for the (111) and (100) faces. The step density is defined to be the number of liquid 
atoms at the interface with at least half their bonds to solid neighbors divided by the 
cross-sectional area of the interface. Figure 2 is an example fitting of a finite difference 
run to Monte Carlo simulations. Spurious peaks found in the average of two Ising model 
simulations is not found in the average of ten runs. Figure 3 gives results for k neq as a 
function of normalized growth velocity in the [111] and [100] directions. At 7 m/s the 
Ising model simulations result in a factor of approximately 2. Aziz et al. found that this 
ratio was about 5 at 1.7 m/s and about 3 at 5 m/s [8]. The ratio ^'neqiiiij/^ne^ioo) found 
by experiment and simulation are in reasonable agreement. Figure 3 (b) plots the data 
for both k neq( in) and k neq ^ 0 o) versus the normalized step velocity. (The normalized step 
velocity is the normalized growth velocity divided by the density of steps.) The results 
for both growth directions fall on a single curve when k neq is plotted as a function of step 
velocity. 

Conclusions 

The Monte Carlo simulations suggest that the determining factor in solute trapping is 
the step velocity at the interface. The Spin-1 Ising model predicts greater incorporation 
of an impurity in Si with k eq = .001 on the (111) face than the (100) face for the same 
growth velocity. The magnitude of the effect is in reasonable agreement with experimental 
results. 
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